### Beginning ####
rm(list = ls())
# analysis study 2
pacman::p_load(multiwayvcov, 
               margins, 
               haven, 
               dplyr, 
               readr, 
               tidyverse, 
               texreg, 
               stargazer, 
               fastDummies, 
               purrr, 
               tidyr, 
               stringr)
# devtools::install_github("ChandlerLutz/starpolishr")
library(starpolishr)

### Import data ####

psych_long2 <- read_dta("./data_clean/psych_clean_long2.dta") %>% select(-RussiaThreat) %>% 
  mutate(RSPEDUC = factor(RSPEDUC),
         INCOM = factor(INCOM),
         oblast = as_factor(oblast),
         narr_topic = as_factor(narr_topic)) %>% 
  rename(StoryTrue = narr_true,
         CRT = zCRT_all,
         CRTUnscaled = crt,
         # AOT = AOT_all,
         AOTUnscaled = AOT_allUnscaled,
         RussianOrientation = russianOrientation,
         RussianOrientationUnscaled = russianOrientation_unscaled,
         RussiaNotAThreat = zthreat_russian,
         Female = RSPSEX,
         Age = RSPAGE,
         Education = RSPEDUC,
         Income = INCOM,
         OblastStrata = strtcode) %>% 
  mutate(Belief = narr_believe,
         AOT = (AOTUnscaled - mean(AOTUnscaled, na.rm = TRUE)) / sd(AOTUnscaled, na.rm = TRUE),
         CRTUnscaled = CRTUnscaled/4,
         AOTUnscaled = (AOTUnscaled - 1)/4)

psych_short2 <- read_csv("./data_clean/psych_clean_short2.csv") %>% 
  rename(CRT = zCRT_all,
         CRTUnscaled = crt.x,
         # AOT = AOT_all,
         AOTUnscaled = AOT_allUnscaled,
         RussianOrientation = russianOrientation,
         RussianOrientationUnscaled = russianOrientation_unscaled,
         RussiaNotAThreat = zthreat_russian,
         RussiaNotAThreatUnscaled = threat_russian,
         Female = RSPSEX,
         Age = RSPAGE,
         Education = RSPEDUC,
         Income = INCOM) %>% 
  mutate(AOT = (AOTUnscaled - mean(AOTUnscaled, na.rm = TRUE)) / sd(AOTUnscaled, na.rm = TRUE),
         CRTUnscaled = CRTUnscaled/4,
         AOTUnscaled = (AOTUnscaled - 1)/4)

# +++++++++++++++++++++++++++++
# Summary Statistics Table ####
# +++++++++++++++++++++++++++++

sum_funs <- list(mean = function(x)mean(x, na.rm = TRUE), median = function(x)median(x, na.rm = TRUE), 
                 min =function(x)min(x, na.rm = TRUE), max = function(x)max(x, na.rm = TRUE),
                 se = function(x)sd(x, na.rm = TRUE)/sqrt(n()), 
                 sd = function(x)sd(x, na.rm = T))

# scaled variables

df.sum <- cbind(
  summarise_all(
    select(psych_short2, CRT, AOT, RussianOrientation, RussiaNotAThreat, Female, Age, Education, Income),
    .funs = sum_funs),
  summarise_all(
    select(psych_long2, Belief, StoryTrue), .funs = sum_funs
  )
)

df.stats.tidy <- 
  df.sum %>% gather(stat, val) %>%
  separate(stat, into = c("var", "stat"), sep = "_") %>%
  spread(stat, val) %>% 
  mutate(var = factor(var, levels = c("Belief", "StoryTrue",
                                      "CRT", "AOT", 
                                      "RussianOrientation",
                                      "RussiaNotAThreat",
                                      "Education", "Income", "Age", "Female"))) %>% 
  select(var, min, mean, median, max, se, sd) %>% # reorder columns
  rename(Variable = var,
         Min. = min,
         Mean = mean,
         Median = median,
         Max. = max, 
         `Std. Error` = se,
         `Std. Dev.` = sd); # rename columns

df.stats.tidy <- df.stats.tidy[order(df.stats.tidy$Variable), ] # reorder manually by factor levels for Variable

t <- xtable::xtable(df.stats.tidy, type = "latex")
print(t, file = "./tables/sumstats_s2.tex", include.rownames = FALSE)

# unscaled variables

df.sum_us <- cbind(
  summarise_all(
    select(psych_short2, CRTUnscaled, AOTUnscaled, RussianOrientationUnscaled, RussiaNotAThreatUnscaled, Female, Age, Education, Income),
    .funs = sum_funs),
  summarise_all(
    select(psych_long2, Belief, StoryTrue), .funs = sum_funs
  )
)

df.stats.tidy_us <- 
  df.sum_us %>% gather(stat, val) %>%
  separate(stat, into = c("var", "stat"), sep = "_") %>%
  spread(stat, val) %>% 
  mutate(var = factor(var, levels = c("Belief", "StoryTrue",
                                      "CRTUnscaled", "AOTUnscaled", 
                                      "RussianOrientationUnscaled",
                                      "RussiaNotAThreatUnscaled",
                                      "Education", "Income", "Age", "Female"))) %>% 
  select(var, min, mean, median, max, se, sd) %>% # reorder columns
  rename(Variable = var,
         Min. = min,
         Mean = mean,
         Median = median,
         Max. = max, 
         `Std. Error` = se,
         `Std. Dev.` = sd); # rename columns

df.stats.tidy_us <- df.stats.tidy_us[order(df.stats.tidy_us$Variable), ] # reorder manually by factor levels for Variable

t_us <- xtable::xtable(df.stats.tidy_us, type = "latex")
print(t_us, file = "./tables/sumstats_s2_unscaled.tex", include.rownames = FALSE)

# ++++++++++++++++++++++++++++
# Zero-Order Correlations ####
# ++++++++++++++++++++++++++++

# full correlations CRT, AOT, and each of the 16 forms of statement

psych_long2 <-
  psych_long2 %>% mutate(
    strat_by_theme = case_when(
      narr_topic == "economic" & narr_strat == "anti-West" ~ "economic|antiwest",
      narr_topic == "economic" & narr_strat == "pro-Russia" ~ "economic|proRussia",
      narr_topic == "economic" & narr_strat == "undermine\nUkraine" ~ "economic|undermineUkraine",
      narr_topic == "economic" & narr_strat == "true story" ~ "economic|truestory",
      narr_topic == "political" & narr_strat == "anti-West" ~ "political|antiwest",
      narr_topic == "political" & narr_strat == "pro-Russia" ~ "political|proRussia",
      narr_topic == "political" & narr_strat == "undermine\nUkraine" ~ "political|undermineUkraine",
      narr_topic == "political" & narr_strat == "true story" ~ "political|truestory",
      narr_topic == "historical" & narr_strat == "anti-West" ~ "historical|antiwest",
      narr_topic == "historical" & narr_strat == "pro-Russia" ~ "historical|proRussia",
      narr_topic == "historical" & narr_strat == "undermine\nUkraine" ~ "historical|undermineUkraine",
      narr_topic == "historical" & narr_strat == "true story" ~ "historical|truestory",
      narr_topic == "military" & narr_strat == "anti-West" ~ "military|antiwest",
      narr_topic == "military" & narr_strat == "pro-Russia" ~ "military|proRussia",
      narr_topic == "military" & narr_strat == "undermine\nUkraine" ~ "military|undermineUkraine",
      narr_topic == "military" & narr_strat == "true story" ~ "military|truestory",
    )
  )

# Full Sample
cortabs <- 
  psych_long2 %>% 
  split(.$strat_by_theme) %>% 
  map(~ cor(select(., CRT, AOT, narr_believe), use = "complete.obs")) %>% 
  data.frame() %>%
  mutate(X = rownames(.)) %>% 
  pivot_longer(-X) %>% 
  separate(name, c("Topic", "Strategy", "Y"), "[.]") %>% 
  filter(value != 1) %>%
  mutate(value = round(value, 3),
         Strategy = factor(Strategy, levels = c("antiwest", 
                                                "proRussia", 
                                                "undermineUkraine", 
                                                "truestory"))) %>% 
  arrange(X, Topic, Strategy) %>%  
  pivot_wider(names_from = Y, values_from = value) %>% 
  split(.$X) %>% 
  lapply(., "[", c(-1))

cortabs$AOT <- cortabs$AOT[, c(-3, -5)]
cortabs$CRT <- cortabs$CRT[, c(-3, -5)]
cortabs$narr_believe <- cortabs$narr_believe[, -4]

stargazer(cortabs$narr_believe %>% 
            mutate(Strategy = as.character(Strategy)), 
          type = "latex", 
          title = "Study 2 Correlations with Belief in Story",
          summary = FALSE, 
          header = FALSE,
          out = "tables/cors1_s2.tex")

# correlations CRT, AOT, Ukraine vs Russia support, having voted, income, education
cordat_s2 <-
  psych_short2 %>%
  dplyr::select(
    CRT,
    AOT,
    RussianOrientation,
    RussiaNotAThreat,
    VOTLOTOM,
    Income,
    Education) %>% 
  mutate(Income = as.numeric(Income),
         Education = as.numeric(Education)) %>% 
  rename("Will Vote" = VOTLOTOM)

cordat_s2 %>%
  cor(use = "complete.obs") %>%
  stargazer(title = "Study 2 Zero-Order Correlations (2)",
            out = "tables/cors2_s2.tex",
            type = "latex",
            header = FALSE)

# +++++++++++++++++++++++++++
# Linear Models - scaled ####
# +++++++++++++++++++++++++++

# part 1

m1 <- list()
m1[[1]] <- lm(narr_believe ~ CRT*StoryTrue, psych_long2)
m1[[2]] <- lm(narr_believe ~ AOT*StoryTrue, psych_long2)

vcov1 <- list()
cl <- makeCluster(4)
vcov1 <- lapply(m1, cluster.vcov, ~ id + narrative, parallel = TRUE)
stopCluster(cl)

# part 2

m2 <- list()
m2[[1]] <- lm(narr_believe ~ CRT*StoryTrue*RussianOrientation, psych_long2)
m2[[2]] <- lm(narr_believe ~ AOT*StoryTrue*RussianOrientation, psych_long2)
m2[[3]] <- lm(narr_believe ~ CRT*StoryTrue*RussiaNotAThreat, psych_long2)
m2[[4]] <- lm(narr_believe ~ AOT*StoryTrue*RussiaNotAThreat, psych_long2)

vcov2 <- list()
cl <- makeCluster(cl)
vcov2 <- lapply(m2, cluster.vcov, ~ id + narrative, parallel = TRUE)
stopCluster(cl)

# part 3

m3 <- list()
m3[[1]] <- lm(narr_believe ~ CRT*StoryTrue*RussianOrientation + 
                Age + Female + factor(Education) + Income + as_factor(OblastStrata), psych_long2)
m3[[2]] <- lm(narr_believe ~ AOT*StoryTrue*RussianOrientation + 
                Age + Female + factor(Education) + Income + as_factor(OblastStrata), psych_long2)
m3[[3]] <- lm(narr_believe ~ CRT*StoryTrue*RussiaNotAThreat + 
                Age + Female + factor(Education) + Income + as_factor(OblastStrata), psych_long2)
m3[[4]] <- lm(narr_believe ~ AOT*StoryTrue*RussiaNotAThreat + 
                Age + Female + factor(Education) + Income + as_factor(OblastStrata), psych_long2)

vcov3 <- list()
cl <- makeCluster(4)
vcov3 <- lapply(m3, cluster.vcov, ~ id + narrative, parallel = TRUE)
stopCluster(cl)

# get SEs
cl.rob.se.1 <- lapply(vcov1, function(x) sqrt(diag(x)))
cl.rob.se.2 <- lapply(vcov2, function(x) sqrt(diag(x)))
cl.rob.se.3 <- lapply(vcov3, function(x) sqrt(diag(x)))

# get F Stats
cl.wald1 <- list()
for(i in 1:length(m1)) {
  cl.wald1[[i]]  <- lmtest::waldtest(m1[[i]], vcov = vcov1[[i]])
}

cl.wald2 <- list()
for(i in 1:length(m2)) {
  cl.wald2[[i]]  <- lmtest::waldtest(m2[[i]], vcov = vcov2[[i]])
}

cl.wald3 <- list()
for(i in 1:length(m3)) {
  cl.wald3[[i]]  <- lmtest::waldtest(m3[[i]], vcov = vcov3[[i]])
}

# get P-Values
p.vals1 <- list()
for(i in 1:length(m1)) {
  p.vals1[[i]] <- lmtest::coeftest(m1[[i]], vcov = vcov1[[i]])[,4]
}

p.vals2 <- list()
for(i in 1:length(m2)) {
  p.vals2[[i]] <- lmtest::coeftest(m2[[i]], vcov = vcov2[[i]])[,4]
}

p.vals3 <- list()
for(i in 1:length(m3)) {
  p.vals3[[i]] <- lmtest::coeftest(m3[[i]], vcov = vcov3[[i]])[,4]
}

# saving main models
main_models_study2 <- list(m1, m2, m3)
save(main_models_study2, file =  "./data_clean/main_clustered_models_study2.Rdata")

# exporting table results
mod.names <- vector("character", sum(length(m1), length(m2), length(m3)))
for(i in 1:sum(length(m1), length(m2), length(m3))) {
  mod.names[i] <- paste0("(", i, ")")
}

tab2 <- 
  stargazer(m1, m2, m3,
            omit = c("Age", "Female", "Education", "Income", "OblastStrata"), 
            # type = "latex",
            column.labels = mod.names,
            type = "latex",
            title = "Study 2 Linear Models",
            dep.var.labels = "Belief",
            add.lines = list(c("Demographic Controls", 
                             rep("No", length(m1)), 
                             rep("No", length(m2)), 
                             rep("Yes", length(m3))
            ),
            c("Oblast Strata Fixed Effects", 
              rep("No", length(m1)), 
              rep("No", length(m2)), 
              rep("Yes", length(m3)))),
            table.layout = "=ldc-ta-s-n",
            se = list(cl.rob.se.1[[1]], cl.rob.se.1[[2]],
                      cl.rob.se.2[[1]], cl.rob.se.2[[2]], cl.rob.se.2[[3]], cl.rob.se.2[[4]],
                      cl.rob.se.3[[1]], cl.rob.se.3[[2]], cl.rob.se.3[[3]], cl.rob.se.3[[4]]),
            p = list(p.vals1[[1]], p.vals1[[2]],
                     p.vals2[[1]], p.vals2[[2]], p.vals2[[3]], p.vals2[[4]],
                     p.vals3[[1]], p.vals3[[2]], p.vals3[[3]], p.vals3[[4]]),
            omit.stat = c("f", "ser"),
            # no.space = T,
            column.sep.width = "-5pt",
            intercept.bottom = FALSE,
            font.size = "tiny",
            order = c(1, 4),
            notes = c("Standard errors are clustered two-ways on Respondent and Narrative. ∗p<0.1;∗∗p<0.05;∗∗∗p<0.01"), 
            notes.align = "l",
            notes.append = FALSE,
            out.header = FALSE,
            header = FALSE
  )  %>% 
  # inserting three stars for each but need to make sure this is true
  star_insert_row(paste0("F Statistic", " ", 
                         "(df = ", 
                         #abs(cl.wald1[[1]][2,]$Df),  ";", " ", 
                         cl.wald1[[1]][2,]$Res.Df, ")", "&", 
                         round(cl.wald1[[1]][2,]$F, 3), "***", "&", 
                         round(cl.wald1[[2]][2,]$F, 3), "***", "&", 
                         # round(cl.wald1[[3]][2,]$F, 3), "***", "&",
                         round(cl.wald2[[1]][2,]$F, 3), "***", "&",
                         round(cl.wald2[[2]][2,]$F, 3), "***", "&",
                         round(cl.wald2[[3]][2,]$F, 3), "***", "&",
                         round(cl.wald2[[4]][2,]$F, 3), "***", '&',
                         round(cl.wald3[[1]][2,]$F, 3), "***", "&",
                         round(cl.wald3[[2]][2,]$F, 3), "***", "&",
                         round(cl.wald3[[3]][2,]$F, 3), "***", "&",
                         round(cl.wald3[[4]][2,]$F, 3), "***",
                         "\\\\"), 
                  insert.after = 73) #%>% cat

write(tab2, "tables/tab2.tex")

# ++++++++++++++++++++++++++++++++++++++++++++++++
# Appendix - Theme and Strategy Fixed Effects ####
# ++++++++++++++++++++++++++++++++++++++++++++++++

psych_long2$strategy <- psych_long2$narr_strat
psych_long2$theme <- psych_long2$narr_topic

# need to create dummies so that we have anti-west baseline
psych_long2$StrategyUndermineUkraine <- ifelse(psych_long2$strategy == "undermine\nUkraine", 1, 0)
psych_long2$StrategyProRussia <- ifelse(psych_long2$strategy == "pro-Russia", 1, 0)

# part 4 (with FEs for theme + strategy)
m4_s2a <- list()
m4_s2a[[1]] <- lm(narr_believe ~ CRT*StoryTrue*RussianOrientation + Age + Female + 
                    as_factor(Education) + Income + as_factor(theme) + 
                    StrategyUndermineUkraine + StrategyProRussia, psych_long2)
m4_s2a[[2]] <- lm(narr_believe ~ AOT*StoryTrue*RussianOrientation + Age + Female + 
                    as_factor(Education) + Income + as_factor(theme) + 
                    StrategyUndermineUkraine + StrategyProRussia, psych_long2)
m4_s2a[[3]] <- lm(narr_believe ~ CRT*StoryTrue*RussiaNotAThreat + Age + Female + 
                    as_factor(Education) + Income + as_factor(theme) + 
                    StrategyUndermineUkraine + StrategyProRussia, psych_long2)
m4_s2a[[4]] <- lm(narr_believe ~ AOT*StoryTrue*RussiaNotAThreat + Age + Female + 
                    as_factor(Education) + Income + as_factor(theme) + 
                    StrategyUndermineUkraine + StrategyProRussia, psych_long2)

vcov4a <- list()
cl <- makeCluster(4)
vcov4a <- lapply(m4_s2a, cluster.vcov, ~ id + narrative, parallel = TRUE)
stopCluster(cl)

# get SEs
cl.rob.se.4a <- lapply(vcov4a, function(x) sqrt(diag(x)))

# get F stats
cl.wald4a <- list()
for(i in 1:length(m4_s2a)) {
  cl.wald4a[[i]]  <- lmtest::waldtest(m4_s2a[[i]], vcov = vcov4a[[i]])
}

# get p-values
p.vals4 <- list()
for(i in 1:length(m4_s2a)) {
  p.vals4[[i]] <- lmtest::coeftest(m4_s2a[[i]], vcov = vcov4a[[i]])[,4]
}

# exporting table results
mod.names <- vector("character", length(m4_s2a))
for(i in 1:length(m4_s2a)) {
  mod.names[i] <- paste0("(", i, ")")
}

reg_m4_s2a <- # strategy baseline = anti-west
  stargazer(m4_s2a,
            #omit = c("Age", "Female", "Education", "Income", "StrategyUndermineUkraine", "StrategyProRussia", "theme"), 
            type = "latex",
            title = "Study 2 Linear Models with Theme and Strategy Controls", 
            column.labels = mod.names,
            # type = "text",
            dep.var.labels = "Belief",
            # add.lines=list(c("Demographic Controls", rep("Yes", length(m4_s2a))),
            #                c("Theme and Strategy Controls", rep("Yes", length(m4_s2a)))), 
            table.layout = "=ldc-ta-s-n",
            se = list( cl.rob.se.4a[[1]], cl.rob.se.4a[[2]], cl.rob.se.4a[[3]], cl.rob.se.4a[[4]]),
            p = list( p.vals4[[1]], p.vals4[[2]], p.vals4[[3]], p.vals4[[4]]),
            omit.stat = c("f", "ser"),
            # no.space = T,
            column.sep.width = "-5pt",
            intercept.bottom = FALSE,
            font.size = "tiny",
            order = c(1, 6),
            notes = "This will be replaced",
            notes.align = "l",
            notes.append = FALSE,
            out.header = FALSE,
            header = FALSE
  ) %>% 
  star_insert_row(paste0("F Statistic", " ", 
                         "(df = ", 
                         #abs(cl.wald1[[1]][2,]$Df),  ";", " ", 
                         cl.wald4a[[1]][2,]$Res.Df, ")", "&", 
                         round(cl.wald4a[[1]][2,]$F, 3), "***", "&",
                         round(cl.wald4a[[2]][2,]$F, 3), "***", "&",
                         round(cl.wald4a[[3]][2,]$F, 3), "***", "&",
                         round(cl.wald4a[[4]][2,]$F, 3), "***",
                         "\\\\"), 
                  # insert.after = 73) # with omitting coefficients
                  insert.after = 116) # no omitting coefficients

note.latex <- "\\multicolumn{5}{l} {\\parbox[t]{11cm}{ \\textit{Notes:} Standard errors are clustered two-ways on Respondent and Narrative. ∗p<0.1;∗∗p<0.05;∗∗∗p<0.0}} \\\\"
reg_m4_s2a[grepl("Note", reg_m4_s2a)] <- note.latex
# cat(reg_m4_s2a, sep = "\n")

write(reg_m4_s2a, file = "tables/reg4_s2a.tex")

rm(psych_long2, psych_short2)

# ++++++++++++++++++++++
# Analysis Unscaled ####
# ++++++++++++++++++++++

psych_long2_us <- read_dta("./data_clean/psych_clean_long2.dta") %>% 
  mutate(RSPEDUC = factor(RSPEDUC),
         INCOM = factor(INCOM),
         oblast = as_factor(oblast),
         narr_topic = as_factor(narr_topic)) %>% 
  rename(StoryTrue = narr_true,
         # CRT = zCRT_all,
         CRT = crt,
         # AOT = AOT_all,
         AOT = AOT_allUnscaled,
         # RussianOrientation = russianOrientation,
         RussianOrientation = russianOrientation_unscaled,
         # RussiaNotAThreat = zthreat_russian,
         # RussiaNotAThreatUnscaled = RussiaNotAThreat,
         RussiaNotAThreat = RussiaThreat,
         Female = RSPSEX,
         Age = RSPAGE,
         Education = RSPEDUC,
         Income = INCOM,
         OblastStrata = strtcode) %>% 
  mutate(Belief = narr_believe,
         OblastStrata = as_factor(OblastStrata),
         CRT = CRT/4,
         AOT = (AOT - 1)/4)


# Linear models

# part 1
m1_us <- list()
m1_us[[1]] <- lm(narr_believe ~ CRT*StoryTrue + OblastStrata, psych_long2_us)
m1_us[[2]] <- lm(narr_believe ~ AOT*StoryTrue + OblastStrata, psych_long2_us)

vcov1_us <- list()
cl <- makeCluster(4)
vcov1_us <- lapply(m1_us, cluster.vcov, ~ id + narrative, parallel = TRUE)
stopCluster(cl)

# part 2
m2_us <- list()
m2_us[[1]] <- lm(narr_believe ~ CRT*StoryTrue*RussianOrientation + OblastStrata, psych_long2_us)
m2_us[[2]] <- lm(narr_believe ~ AOT*StoryTrue*RussianOrientation + OblastStrata, psych_long2_us)
m2_us[[3]] <- lm(narr_believe ~ CRT*StoryTrue*(RussiaNotAThreat) + OblastStrata, psych_long2_us)
m2_us[[4]] <- lm(narr_believe ~ AOT*StoryTrue*(RussiaNotAThreat) + OblastStrata, psych_long2_us)

vcov2_us <- list()
cl <- makeCluster(4)
vcov2_us <- lapply(m2_us, cluster.vcov, ~ id + narrative, parallel = TRUE)
stopCluster(cl)

# saving main models as .rds
write_rds(m2_us, file = "./data_clean/study2_main_models_unscaled.rds")

# part 3
m3_us <- list()
m3_us[[1]] <- lm(narr_believe ~ CRT*StoryTrue*RussianOrientation + 
                Age + Female + factor(Education) + Income + OblastStrata, psych_long2_us)
m3_us[[2]] <- lm(narr_believe ~ AOT*StoryTrue*RussianOrientation + 
                Age + Female + factor(Education) + Income + OblastStrata, psych_long2_us)
m3_us[[3]] <- lm(narr_believe ~ CRT*StoryTrue*RussiaNotAThreat + 
                Age + Female + factor(Education) + Income + OblastStrata, psych_long2_us)
m3_us[[4]] <- lm(narr_believe ~ AOT*StoryTrue*RussiaNotAThreat + 
                Age + Female + factor(Education) + Income + OblastStrata, psych_long2_us)

vcov3_us <- list()
cl <- makeCluster(4)
vcov3_us <- lapply(m3_us, cluster.vcov, ~ id + narrative)
stopCluster(cl)

# get SEs
cl_rob_se_1_us <- lapply(vcov1_us, function(x) sqrt(diag(x)))
cl_rob_se_2_us <- lapply(vcov2_us, function(x) sqrt(diag(x)))
cl_rob_se_3_us <- lapply(vcov3_us, function(x) sqrt(diag(x)))

# get F Stats
cl_wald1_us <- list()
for(i in 1:length(m1_us)) {
  cl_wald1_us[[i]]  <- lmtest::waldtest(m1_us[[i]], vcov = vcov1_us[[i]])
}

cl_wald2_us <- list()
for(i in 1:length(m2_us)) {
  cl_wald2_us[[i]]  <- lmtest::waldtest(m2_us[[i]], vcov = vcov2_us[[i]])
}

cl_wald3_us <- list()
for(i in 1:length(m3_us)) {
  cl_wald3_us[[i]]  <- lmtest::waldtest(m3_us[[i]], vcov = vcov3_us[[i]])
}

# get P-Values
pvals1_us <- list()
for(i in 1:length(m1_us)) {
  pvals1_us[[i]] <- lmtest::coeftest(m1_us[[i]], vcov = vcov1_us[[i]])[,4]
}

pvals2_us <- list()
for(i in 1:length(m2_us)) {
  pvals2_us[[i]] <- lmtest::coeftest(m2_us[[i]], vcov = vcov2_us[[i]])[,4]
}

pvals3_us <- list()
for(i in 1:length(m3_us)) {
  pvals3_us[[i]] <- lmtest::coeftest(m3_us[[i]], vcov = vcov3_us[[i]])[,4]
}

# saving main models
main_models_study2_us <- list(m1_us, m2_us, m3_us)
save(main_models_study2_us,file =  "./data_clean/main_clustered_models_study2_unscaled.Rdata")

# exporting table results
mod_names_us <- vector("character", sum(length(m1_us), length(m2_us), length(m3_us)))
for(i in 1:sum(length(m1_us), length(m2_us), length(m3_us))) {
  mod_names_us[i] <- paste0("(", i, ")")
}

tab2 <- 
  stargazer(m1_us, m2_us, m3_us,
            omit = c("Age", "Female", "Education", "Income", "OblastStrata"), 
            # type = "latex",
            column.labels = mod_names_us,
            type = "latex",
            title = "Study 2 Linear Models",
            dep.var.labels = "Belief",
            add.lines = list(c("Demographic Controls", 
                               rep("Yes", length(m1_us)), 
                               rep("Yes", length(m2_us)), 
                               rep("Yes", length(m3_us))
            ),
            c("Oblast Strata Fixed Effects", 
              rep("No", length(m1_us)), 
              rep("No", length(m2_us)), 
              rep("Yes", length(m3_us)))),
            table.layout = "=ldc-ta-s-n",
            se = list(cl_rob_se_1_us[[1]], cl_rob_se_1_us[[2]],
                      cl_rob_se_2_us[[1]], cl_rob_se_2_us[[2]], cl_rob_se_2_us[[3]], cl_rob_se_2_us[[4]],
                      cl_rob_se_3_us[[1]], cl_rob_se_3_us[[2]], cl_rob_se_3_us[[3]], cl_rob_se_3_us[[4]]),
            p = list(pvals1_us[[1]], pvals1_us[[2]],
                     pvals2_us[[1]], pvals2_us[[2]], pvals2_us[[3]], pvals2_us[[4]],
                     pvals3_us[[1]], pvals3_us[[2]], pvals3_us[[3]], pvals3_us[[4]]),
            omit.stat = c("f", "ser"),
            # no.space = T,
            column.sep.width = "-5pt",
            intercept.bottom = FALSE,
            font.size = "tiny",
            order = c(1, 4),
            notes = c("Standard errors are clustered two-ways on Respondent and Narrative. ∗p<0.1;∗∗p<0.05;∗∗∗p<0.01"), 
            notes.align = "l",
            notes.append = FALSE,
            out.header = FALSE,
            header = FALSE
  )  %>% 
  # inserting three stars for each but need to make sure this is true
  star_insert_row(paste0("F Statistic", " ", 
                         "(df = ", 
                         #abs(cl_wald1_us[[1]][2,]$Df),  ";", " ", 
                         cl_wald1_us[[1]][2,]$Res.Df, ")", "&", 
                         round(cl_wald1_us[[1]][2,]$F, 3), "***", "&", 
                         round(cl_wald1_us[[2]][2,]$F, 3), "***", "&", 
                         # round(cl_wald1_us[[3]][2,]$F, 3), "***", "&",
                         round(cl_wald2_us[[1]][2,]$F, 3), "***", "&",
                         round(cl_wald2_us[[2]][2,]$F, 3), "***", "&",
                         round(cl_wald2_us[[3]][2,]$F, 3), "***", "&",
                         round(cl_wald2_us[[4]][2,]$F, 3), "***", '&',
                         round(cl_wald3_us[[1]][2,]$F, 3), "***", "&",
                         round(cl_wald3_us[[2]][2,]$F, 3), "***", "&",
                         round(cl_wald3_us[[3]][2,]$F, 3), "***", "&",
                         round(cl_wald3_us[[4]][2,]$F, 3), "***",
                         "\\\\"), 
                  insert.after = 73) #%>% cat

write(tab2, "tables/tab2_unscaled.tex")

tab4 <- 
  stargazer(m1_us, m2_us, m3_us,
            omit = c("Age", "Female", "Education", "Income", "OblastStrata"), 
            # type = "latex",
            column.labels = mod_names_us,
            type = "html",
            title = "Study 2 Linear Models",
            dep.var.labels = "Belief",
            add.lines = list(c("Demographic Controls", 
                               rep("No", length(m1_us)), 
                               rep("No", length(m2_us)), 
                               rep("Yes", length(m3_us))
            ),
            c("Oblast Strata Fixed Effects", 
              rep("Yes", length(m1_us)), 
              rep("Yes", length(m2_us)), 
              rep("Yes", length(m3_us)))),
            table.layout = "=ldc-ta-s-n",
            se = list(cl_rob_se_1_us[[1]], cl_rob_se_1_us[[2]],
                      cl_rob_se_2_us[[1]], cl_rob_se_2_us[[2]], cl_rob_se_2_us[[3]], cl_rob_se_2_us[[4]],
                      cl_rob_se_3_us[[1]], cl_rob_se_3_us[[2]], cl_rob_se_3_us[[3]], cl_rob_se_3_us[[4]]),
            p = list(pvals1_us[[1]], pvals1_us[[2]],
                     pvals2_us[[1]], pvals2_us[[2]], pvals2_us[[3]], pvals2_us[[4]],
                     pvals3_us[[1]], pvals3_us[[2]], pvals3_us[[3]], pvals3_us[[4]]),
            omit.stat = c("f", "ser"),
            # no.space = T,
            column.sep.width = "-5pt",
            intercept.bottom = FALSE,
            font.size = "tiny",
            order = c(1, 4),
            notes = c("Standard errors are clustered two-ways on Respondent and Narrative. ∗p<0.1;∗∗p<0.05;∗∗∗p<0.01"), 
            notes.align = "l",
            notes.append = FALSE,
            out.header = FALSE,
            header = FALSE
  )  %>% 
  # inserting three stars for each but need to make sure this is true
  star_insert_row(paste0("<tr><td style=\"text-align:left\">", 
                         "F Statistic", " ", 
                         "(df = ", 
                         #abs(cl_wald1_us[[1]][2,]$Df),  ";", " ", 
                         cl_wald1_us[[1]][2,]$Res.Df, ")", "</td><td>", 
                         round(cl_wald1_us[[1]][2,]$F, 3), "***", "</td><td>", 
                         round(cl_wald1_us[[2]][2,]$F, 3), "***", "</td><td>", 
                         # round(cl_wald1_us[[3]][2,]$F, 3), "***", "</td><td>",
                         round(cl_wald2_us[[1]][2,]$F, 3), "***", "</td><td>",
                         round(cl_wald2_us[[2]][2,]$F, 3), "***", "</td><td>",
                         round(cl_wald2_us[[3]][2,]$F, 3), "***", "</td><td>",
                         round(cl_wald2_us[[4]][2,]$F, 3), "***", '</td><td>',
                         round(cl_wald3_us[[1]][2,]$F, 3), "***", "</td><td>",
                         round(cl_wald3_us[[2]][2,]$F, 3), "***", "</td><td>",
                         round(cl_wald3_us[[3]][2,]$F, 3), "***", "</td><td>",
                         round(cl_wald3_us[[4]][2,]$F, 3), "***", "</td></tr>"),
                         #"\\\\"), 
                  insert.after = 73) #%>% cat

write(tab4, file = "./tables/tab2_unscaled.html")
